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Abstract:  COMSOL  Multiphysics  was  used  to 
solve  a  phonon  Boltzmann  transport  equation 
(BTE)  for  nanoscale  heat  transport  problems. 
One-dimensional  steady-state  and  transient  BTE 
problems  were  successfully  solved  based  on 
finite  element  and  discrete  ordinate  methods  for 
spatial  and  angular  discretizations,  respectively, 
by  utilizing  the  built-in  feature  of  the  COMSOL, 
Coefficient  Form  of  PDF.  A  sensitivity  study 
was  conducted  with  various  discretization 
refinements  for  different  values  of  the  Knudsen 
number,  which  is  a  measure  of  the  nanoscale 
regime.  It  was  found  that  sufficient  refinement 
for  angular  discretization  is  critical  in  obtaining 
accurate  solutions  of  the  BTE. 

Keywords:  Nanoscale  heat  transport,  phonon 
Boltzmann  transport  equation,  finite  element 
method,  discrete  ordinates  method. 

1.  Introduction 

For  the  last  two  centuries,  many  heat 
conduction  problems  have  been  modeled 
successfully  by  a  Fourier  diffusive  equation 
(FDE).  The  equation  can  be  derived  using  a 
conservation  law  of  energy  and  Fourier’s  linear 
approximation  of  heat  flux  using  a  temperature 
gradient.  The  FDE  is  a  parabolic  equation 
reflecting  a  diffusive  nature  of  heat  transport. 
An  underlying  assumption  is  that  the  heat  is 
effectively  transferred  between  localized  regions 
through  sufficient  scattering  events  of  phonons 
within  a  medium.  Therefore,  the  FDE  does  not 
hold  when  the  number  of  scatterings  is  negligible, 
which  could  happen  when  a  mean  free  path  of 
phonon  is  similar  or  a  larger  order  of  magnitude 
than  a  domain  size  of  interest,  significant 
boundary  scattering  occurs  at  material  interfaces 
causing  thermal  resistance,  etc.  Another 
problem  of  the  FDE  is  that  it  admits  an  infinite 
speed  of  the  heat  transport,  which  is 
contradictory  to  the  theory  of  relativity. 
Therefore,  the  FDE  is  inappropriate  for  the  heat 
transfer  problems  at  small  time  and  spatial  scales. 


To  resolve  the  issue  of  the  infinite  speed  of 
heat  carrier,  hyperbolic  equationns,  such  as  a 
Cattaneo  equation  (Telegraph  equation)  [1]  or  a 
relativistic  heat  equation  [2],  are  often  used  to 
reflect  a  wave  nature  of  the  heat  transport  for 
cases  of  small  speed  of  heat  propagation  (e.g., 
low  temperature,  low  conductive  medium, 
thermal  insulator,  phase  transition,  etc.)  and/or 
high  speed  of  heat  flux  (e.g.,  pico-  or  femto¬ 
second  pulsed-laser  heating.)  Despite  some 
issues  [1],  these  hyperbolic  equations  can  be 
used  to  consider  the  finite  speed  of  phonons  for  a 
short  time  scale.  However,  they  still  cannot  be 
used  for  small  spatial  scale. 

For  the  nanoscale  heat  transfer  analysis,  one 
needs  equations  and  methods  for  small-scale 
simulation  in  terms  of  both  time  and  space.  A 
molecular  dynamics  (MD)  simulation  can  be  a 
useful  and  accurate  method  in  this  regard. 
However,  the  MD  simulation  is  usually 
computationally  expensive,  so  it  is  suitable  for 
systems  having  a  few  atomic  layers  or  several 
thousand  atoms,  but  not  suitable  for  device-level 
thermal  analysis.  Alternatively,  a  phonon 
Boltzmann  transport  equation  (BTE),  a  first- 
order  partial  differential  equation  for  a  phonon 
distribution  function,  has  been  used  for  this 
purpose  [3-5].  The  distribution  function  is  a 
scalar  quantity  in  the  six-dimensional  phase 
space  (three  space  coordinates  and  three  wave- 
vector  coordinates).  The  phonon  BTE  is  also 
called  an  equation  of  phonon  radiative  transfer 
(EPRT)  when  the  phonon  distribution  function  is 
replaced  with  a  phonon  intensity  function  [3]. 
The  BTE  is  known  to  be  difficult  to  solve,  and 
thus  often  simplified  with  a  relaxation  time 
approximation.  The  BTE  can  predict  a  ballistic 
nature  of  heat  transfer  under  an  assumption  that 
particle-like  behaviors  of  phonon  are  much  more 
significant  than  its  wave-like  behaviors,  which 
make  the  BTE  valid  for  structures  larger  than  the 
wavelength  of  phonons.  The  BTE  can  be  solved 
analytically  for  simple  geometries  [4,  6,  7],  or 
numerically  for  complex  geometries  by  using 
either  deterministic  methods  (e.g.,  discrete 
ordinates  method,  spherical  harmonics  method, 


Report  Documentation  Page 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 


1.  REPORT  DATE 

OCT  2009 


2.  REPORT  TYPE 


4.  TITLE  AND  SUBTITLE 

Nanoscale  Heat  Transfer  using  Phonon  Boltzmann  Transport  Equation 


6.  AUTHOR(S) 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory ,300  College  Park, Dayton, OH, 45469-0060 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


3.  DATES  COVERED 

00-00-2009  to  00-00-2009 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

5d.  PROIECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

Presented  at  the  2009  COMSOL  Conference  Boston,  8-10  Oct,  Boston,  MA 

14.  ABSTRACT 

COMSOL  Multiphysics  was  used  to  solve  a  phonon  Boltzmann  transport  equation  (BTE)  for  nanoscale 
heat  transport  problems.  One-dimensional  steady-state  and  transient  BTE  problems  were  successfully 
solved  based  on  finite  element  and  discrete  ordinate  methods  for  spatial  and  angular  discretizations, 
respectively  by  utilizing  the  built-in  feature  of  the  COMSOL  Coefficient  Form  of  PDE.  A  sensitivity  study 
was  conducted  with  various  discretization  refinements  for  different  values  of  the  Knudsen  number,  which 
is  a  measure  of  the  nanoscale  regime.  It  was  found  that  sufficient  refinement  for  angular  discretization  is 
critical  in  obtaining  accurate  solutions  of  the  BTE. 

15.  SUBIECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

19a.  NAME  OF 

ABSTRACT 

OF  PAGES 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

23 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


finite  volume  method  [8],  etc.)  or  statistical 
methods  (e.g.,  Monte  Carlo  simulation  [9,  10]). 
In  general,  solving  the  BTE  is  much  more 
efficient  than  the  MD,  and  the  predictions  agree 
well  with  experimental  data  [10,  11]. 

In  the  present  study,  we  will  solve  the  BTE 
with  the  relaxation  time  approximation  for  one¬ 
dimensional  heat  transport  problems  using  the 
COMSOL  Multiphysics.  We  will  combine  the 
finite  element  (FE)  capability  of  the  COMSOL 
with  the  DOM  to  solve  both  steady-state  and 
transient  heat  transfer  problems.  The  numerical 
results  will  be  compared  with  analytical  ones  for 
the  steady-state  problem.  A  sensitivity  study 
will  be  conducted  for  various  scale  measures  to 
show  the  difference  between  the  BTE  solutions 
with  the  conventional  FDE  ones.  Furthermore, 
we  will  determine  a  suitable  discretization 
scheme  for  spatial  and  angular  spaces  in  order  to 
obtain  reliable  transient  solutions  in  using  the 
COMSOL. 


2.  Use  of  COMSOL  Multiphysics 


The  BTE  equation  under  a  relaxation  time 
approximation  is  written  as 


3/ 


^r+v-v/  =  K- 


dt 


df 


dt 


fo-f 


(1) 


where  f  is  the  carrier  distribution  function,  f ’ 


the  equilibrium  distribution,  v  the  carrier  group 
velocity  vector,  and  r„  the  relaxation  time.  The 

EPRT  in  terms  of  the  phonon  intensity  is  written 
as 


dlw 

dt 


-  v  ■  VI  - 


where  the  phonon  intensity  is 
I  w(t  ,v  ,r)  =  \v\haf  {t  ,v  ,r)D(co)  !  4jt  , 


(2) 

expressed  as 
where  h  is 


the  reduced  Planck  constant,  co  the  phonon 
angular  frequency,  and  D(co)  the  density  of  state 


per  unit  volume.  An  equilibrium  phonon 
intensity  in  Eq.  (2)  is  written  as 


Im 


4k 


\imdO. 


(3) 


For  one-dimensional  problems,  Eq.  (2)  becomes 


dla 

dt 


-  vcos 0 


dlw 

dx 


(4) 


where  0  is  a  polar  angle  between  a  phonon 
propagation  direction  and  the  global  heat  transfer 
direction  ( x ). 


A  solution  process  of  the  1-D  BTE  in  Eq.  (4) 
along  with  Eq.  (3)  can  be  summarized  as 
follows:  With  an  initial  value  of  /  ,  Eq.  (4)  is 

solved  for  the  phonon  intensity  field,  I  (t,x) ,  for 
a  particular  phonon  propagation  angle,  0 .  After 
obtaining  the  solutions  for  all  possible  solid 
angles,  4x ,  I m  can  be  calculated  using  Eq.  (3). 
With  this  updated  value  of  /  ,  the  solution 


process  repeats  until  the  convergence  of  / '  and 
/  x .  Therefore,  the  BTE  and  thus  the  EPRT  are 


highly  nonlinear  equations  to  solve.  Once  the 
phonon  intensity  is  solved,  heat  flux  and  an 
equivalent  temperature  can  be  obtained  by 
integrating  the  phonon  intensity  over  the  solid 
angles,  4 n ,  such  that 


cos  0  dQ. 

and 

(5) 

4n  Im 

(6) 

C  V  | 

cM  , 

where  c  is  the  volumetric  specific  heat. 

One  of  the  most  widely  used  methods  for 
solving  the  BTE  numerically  is  the  discrete 
ordinates  method  (DOM).  The  DOM  is  a  tool  to 
transform  the  equation  of  radiative  transfer  into  a 
set  of  simultaneous  partial  differential  equations. 
This  is  based  on  a  discrete  representation  of 
directional  variation  of  intensity.  A  solution  to 
the  transport  problem  can  be  found  by  solving 
the  equation  of  radiative  transfer  for  a  set  of 
discrete  directions  spanning  the  entire  solid  angle. 
The  integrals  over  the  solid  angle  can  be 
approximated  by  numerical  integration  using 
Gauss-Legender  quadratures. 

The  BTEs  in  Eqs.  (3)-(4)  were  solved  by  the 
DOM  using  a  built-in  feature  of  Coefficient 
Form  PDEs  in  COMSOL  Multiphysics.  The 
spatial  domain  (  x  )  was  discretized  into  n 
elements  as  a  normal  FE  mesh  refinement,  while 
the  angular  domain  ( Q ),  referred  to  as  ordinates, 
was  discretized  in  the  directional  coordinates  by 
dividing  the  angular  space  into  a  finite  number, 
m  .  The  angular  domain  was  discretized  at 
Gaussian  quadrature  points  within 
- 1  <  n  =  cos  9  <  1  •  For  each  polar  angle  ( 0m ),  we 


solved  the  system  of  nonlinear  equations  for  the 
phonon  intensity,  f  t(t,x,/um ),  and  then  updated 

the  equilibrium  phonon  intensity,  I :  ft,x),  using 


/  for  all  angles.  Eqs.  (3)  and  (5)  can  be  written 
in  discrete  forms  as 

=  and  (7) 

^  m 

q{t,x)=27rY,Iait,x,/im)^mwm 

m 

respectively.  The  weight  satisfies  V  w  =2- 

The  solver  used  for  the  steady-state  and 
transient  analyses  was  a  built-in  direct  solver, 
UMFPACK.  Default  values  were  used  for  the 
solver  parameters  except  that  strict  time  steps 
were  taken  by  the  solver  and  a  maximum  BDF 
order  was  set  to  1.  When  the  relative  error  of  the 
calculated  value  of  the  equivalent  equilibrium 
intensity  between  two  iteration  steps  was  less 
than  a  present  tolerance,  we  considered  the 
problem  was  converged  and  then  calculated  the 
equivalent  temperature  and  heat  flux  using  Eqs. 
(6)  and  (8),  respectively. 

5.  Results  and  Discussion 

We  have  solved  the  BTE  equations  with  the 
COMSOE  under  both  steady-state  and  transient 
states.  Spatial  and  angular  domains  were 
discretized  into  n  elements  and  m  quadrature 
points,  respectively.  A  quadratic  Lagrange 
element  was  used  for  the  spatial  FE  mesh.  The 
total  degree  of  freedom  for  the  1-D  problem 
becomes  2m(n  +  \/2)  .  For  the  steady-state 
problems,  it  was  observed  that  the  well- 
converged  solutions  were  obtained  with  a 
relatively  coarse  mesh  (15  elements)  while 
using  16  quadrature  points.  The  sensitivity  of 
mesh  and  quadrature  points  for  the  transient 
analysis  will  be  stated  later  in  this  section. 

Figure  1  shows  steady-state  temperature 
distributions  via  nondimensional  emissive  power 
(e*  )  along  the  1-D  domain  for  various  Knudsen 

numbers  (  Kn  )•  For  comparison  purpose,  an 
analytical  solution  was  also  obtained  by  using  an 
analogy  to  the  radiation  heat  transfer  [5]  by 
solving  an  integral  equation, 

2 e*(7)  =  E2(n)  +  £  e„ tn')  E, (\q  -  q'\) d q'  (9) 

where  q  =  xj L  is  the  nondimensional  coordinate, 
where  L  is  the  length  of  the  domain, 
En (x)  =  J’ n"~-  exp(-  x/ju)dju  the  exponential 
integral,  £  the  optical  thickness,  which  is  an 


inverse  of  the  Knudsen  number,  Kn  =  A/L  , 
where  A  is  the  average  mean  free  path  of 
phonons,  and  eo*  (q)  =  [eo  (q)  -  J~2  ] /[J+ql  -  J~2  ]  is 

the  emissive  power  nondimensionalized  by  a 
difference  of  boundary  heat  fluxes  ( j~2  and  J+) 

at  both  ends  of  the  1-D  domain.  In  Figure  1,  the 
analytic  and  numerical  solutions  for  various  Kn 
(=  0.1,  1,  10,  100)  were  drawn  with  markers  and 
lines,  respectively,  which  shows  an  excellent 
agreement  of  the  solutions  with  each  other. 
While  the  conventional  Fourier  solution  would 
yield  the  temperature  gradient  varying  from  1  at 
q  =  0  to  0  at  //  =  1 ,  the  solution  from  the  BTE 
yields  a  reduced  temperature  gradient.  The 
larger  the  Knudsen  number,  the  smaller  the 
temperature  gradient,  which  is  attributed  from 
more  ballistic  phonon  transport  as  compared  to 
the  diffusive  one. 


Nondimensional  coordinate,  rj 


Figure  1.  Steady-state  nondimensionalized  emissive 
power  distribution  for  various  Knudsen  numbers. 

We  studied  the  effects  of  refinements  in 
terms  of  spatial  FE  meshes  and  angular 
quadrature  points  on  the  time-dependent 
transient  BTE  solutions.  The  transient  time  ( t ) 
was  nondimensionalized  by  the  relaxation  time 
of  phonons  (t0),  so  that  {  =  t/r„  •  Figure  2 

shows  nondimensionalized  temperature 
distributions  at  t*  =0.1  along  the  1-D  domain  for 
Kn  =10  ■  We  used  three  different  FE  meshes  (15, 
60  and  120  elements)  and  a  fixed  angular 
division  with  1 6  Gaussian  points  ( „sp  =  16  )■  The 

three  FE  meshes  yield  nearly  identical 
temperature  solutions  with  a  decreasing  trend 
from  the  hot  side  ( q  =  o )  to  the  cold  side  ( q  =  l ). 


Smooth  solutions  were  obtained  except  at  a 
region  near  the  hot  boundary.  An  inset 
rectangular  area  in  Figure  2  within  0  <  t)  <  0.3 
was  zoomed  in  and  replotted  in  Figure  3,  which 
clearly  shows  significant  wiggles  near  the  hot 
boundary.  Therefore,  the  FE  mesh  refinement 
does  not  help  smoothing  out  the  wiggles  at  all. 
The  wiggling  of  the  solution,  called  a  ray  effect, 
is  known  to  be  attributed  to  insufficient 
refinement  in  the  angular  domain,  and  thus 
cannot  be  smoothed  out  with  the  FE  mesh 
refinement  in  the  spatial  domain. 


Nondime  ns  ional  coordinate,  r\ 


Figure  2.  Transient  nondimensional  temperature 
distribution  predicted  with  15,  60  and  120  finite 
elements  and  16  Gaussian  points. 


Figure  3.  A  magnified  view  of  Figure  2  within 
0  <  //  <  0.3  ■ 

To  eliminate  the  ray  effect,  we  made  further 
refinement  in  the  angular  domain  by  using  more 
quadrature  points.  Figure  4  shows 
nondimensionalized  temperature  distributions  at 


t‘  =0.1  along  the  1-D  domain  for  Kn  =  10 
predicted  with  six  different  Gaussian  quadrature 
points  („gp  =  4,  8, 16,  32,  64, 128)-  We  utilized  a 

Matlab-scripting  feature  of  COMSOL  to 
implement  the  large  n„p  in  building  the  system 

of  equations.  A  fixed  spatial  division  with  240 
FE  elements  was  used  in  these  calculations.  An 
inset  rectangular  area  within  0  <  rj  <  0.3  was 
magnified  and  replotted  in  Figure  5.  Although  a 
highly  refined  FE  mesh  (240  elements)  was  used, 
the  coarse  angular  divisions  alleviate  the 
wiggling  of  the  solution  (e.g.,  n„p  =  4  )•  The 

wiggles  near  the  hot  end  gradually  decrease  with 
the  increase  of  ,  and  thus  the  ray  effect. 

Therefore,  we  found  that  the  spatial  and  angular 
refinements  are  independent  of  each  other,  and 
thus  do  not  recommend  to  use  the  highly  refined 
spatial  mesh  with  the  coarse  angular  mesh. 


Figure  4.  Transient  nondimensional  temperature 
distribution  predicted  with  240  finite  elements  and  six 
Gaussian  quadrature  points  (  ngp  )  for  angular 

discretization. 


Figure  5.  A  magnified  view  of  Figure  5  within 
0  <  /;  <  0.3  ■ 

Figure  6  and  Figure  7  show  the  temperature 
and  heat  flux  distributions,  respectively,  at 
various  time  scales  with  three  different  Knudsen 
numbers.  The  result  successfully  reproduced  the 
one  reported  earlier  by  Chen  [4],  but  this  time 
using  the  FE  method  using  the  COMSOL.  The 
figures  show  the  heat  transfer  from  the  hot  side 
to  the  cold  side  with  the  finite  speed  of  phonon 
as  the  short  term  (small  f )  solutions  indicate, 
which  cannot  be  predicted  by  the  FDE.  While 
constant  temperatures  were  applied  as  boundary 
conditions  at  both  sides,  we  can  observe  the 
temperature  jump  at  these  boundaries.  The 
higher  the  Kn ,  the  larger  the  jump.  The  reason 
for  the  temperature  jump  at  these  emissive 
boundaries  was  well  explained  in  earlier  work 
[7].  For  a  small  Knudsen  number  at  Kn  =0.1, 
the  temperature  solution  approaches  that  of  the 
steady-state  FDE  as  the  time  increases,  while  for 


a  large  Knudsen  number  at  Kn  =  10  ,  the 
temperature  solution  approaches  a  nearly 
horizontal  line,  which  indicates  the  ballistic  heat 
transfer  rather  than  the  diffusive  one,  and  results 
in  a  small  temperature  gradient  and  thus  a  large 
thermal  conductivity. 

6.  Conclusions 

We  have  solved  the  phonon  Boltzmann 
transport  equation  for  the  nanoscale  heat  transfer 
problems  with  the  COMSOL  Multiphysics. 
One-dimensional  steady-state  and  transient 
problems  were  successfully  solved  using  FEM 
and  DOM  for  spatial  and  angular  discretizations, 
respectively,  by  utilizing  the  built-in  feature  of 
the  COMSOL,  Coefficient  Form  ofPDE. 

A  sensitivity  study  was  conducted  with 
various  discretization  refinements  for  different 
values  of  the  Knudsen  number,  which  is  a 
measure  of  the  nanoscale  regime.  Significant  ray 
effects  were  found  with  coarse  angular 
discretizations.  However,  false  scattering  due  to 
coarse  spatial  discretization  was  not  severe  with 
the  COMSOL  solution.  It  was  found  that  the 
spatial  and  angular  refinements  were 
independent  of  each  other,  and  it  was  not 
recommended  to  use  the  highly  refined  spatial 
mesh  with  the  coarse  angular  mesh  in  solving  the 
BTE  with  the  COMSOL.  The  present  success 
with  the  1  -D  simulation  encourages  us  to  apply  it 
for  multidimensional  cases. 


Figure  6.  Temperature  distributions  at  different  time  scales,  (a)  Kn  =  10  ,  (b)  Kn  =  1  and  (c)  Kn  =  0.1. 


Figure  7.  Heat  flux  distributions  at  different  time  scales,  (a)  Kn  =  10  ,  (b)  Kn  =  1  and  (c)  Kn  =  0.1  ■ 
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•  Background  information. 

•  Description  of  phonon  Boltzmann  transport  equation  (BTE). 
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Summary  and  conclusions. 


Fourier  Equation  (FE) 


For  last  two  centuries,  heat  conduction  has  been  modeled  by  Fourier  Eq  (FE). 

dT  „ 

—  Conservation  of  energy:  pc  —  =  -V  •  q 

dt  Diffusive 

—  Fourier’s  linear  approximation  of  heat  flux:  =  -kVT 

1  dT 


> 


a  dt 


=  VZT 


Parabolic  equation  — >  Diffusive  nature  of  heat  transport. 

Heat  is  effectively  transferred  between  localized  regions  through  sufficient 
scattering  events  of  phonons  within  medium. 

Does  not  hold  when  number  of  scattering  is  negligible. 

—  e.g.,  mean  free  path  ~  device  size  (chip-package  level). 

—  Boundary  scattering  at  interfaces  causing  thermal  resistance. 

Admits  infinite  speed  of  heat  transport  — >  Contradict  with  theory  of  relativity. 


^  Fourier  Equation  cannot  be  used  for  small  time  and  spatial  scales. 
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Hyperbolic  Heat  Conduction  Equation  (HHCE) 


Resolve  the  issue  of  the  Fourier  equation  with  the  infinite  speed  of  heat  carrier. 

> 


1  d2T 

1  dT 

I 

C 2  dt2 

1 

a  dt 

=  VT 


—  Definition  of  heat  flux: 


(C  =a/T0) 


f  q  =  -kV  T  ,  (t0:  relaxation  time) 


•  Hyperbolic  equation  — >  Wave  nature  of  heat  transport. 

•  Called  as  Cattaneo  equation  or  Telegraph  equation. 

•  Finite  speed  of  heat  carriers. 

•  Ad  hoc  approximation  of  heat  flux  definition. 

•  Violates  2nd  law  of  thermodynamics. 

—  If  heat  source  varies  faster  than  speed  of  sound,  heat  would  appear  to  be  moving  from 
cold  to  hot. 


'r  HHCE:  could  be  used  for  short  time  scale,  but  not  for  short  spatial  scale. 


Small  Scale  Heat  Transport  (Time  &  Space) 


•  Fourier  Equation  cannot  be  used  for  small  time  and  spatial  scales. 

•  HHCE:  could  be  used  for  short  time  scale,  but  not  for  short  spatial  scale. 

•  Needs  equations  and  methods  for  small  scale  simulation  in  terms  of  both 
time  and  space. 

—  Molecular  dynamics  simulation. 

•  Accurate  method. 

•  Computationally  expensive. 

•  Suitable  for  systems  having  a  few  atomic  layers  or  several  thousands  of  atoms. 

•  Not  suitable  for  device-level  thermal  analysis. 

—  Boltzmann  Transport  Equation  (BTE). 

—  Ballistic-Diffusive  Equation  (BDE). 

•  Similar  to  Cattaneo  Eq.  (HHCE)  with  a  source  term. 

•  Derived  from  BTE. 

•  Good  approximation  of  BTE  without  internal  heat  source,  disturbance,  etc. 


Boltzmann  Transport  Equation  (BTE) 


•  BTE:  also  called  as  equation  of  phonon  radiative  transfer  (EPRT). 

•  Equation  for  phonon  distribution  function: 


fo~f 


Ballistic 


L«  A  j2 


V  J  scat 


T 


L 


•  Can  predict  ballistic  nature  of  heat  transfer. 

•  Neglects  wave-like  behaviors  of  phonon. 


T. 


—  Valid  for  structures  larger  than  wavelength  of  phonons  (~  1  nm  @  RT). 

•  Solution  methods: 

—  Deterministic:  discrete  ordinates  method,  spherical  harmonics  method. 
—  Statistical:  Monte  Carlo. 

•  Much  more  efficient  than  MD. 

•  Agrees  well  with  experimental  data. 


Details  of  Boltzmann  Transport  Equation  (BTE) 


Phonon  intensity:  Im(t,v,r)  =  \v\ha)f(t,v,r)D(w)/4x 


r\T  J  _ J 

BTE  becomes  EPRT:  — — +  v-V7  =  — — — .  = 


CO 


For  1-D, 


dl 


CO 


dt 


+  vcos# 


dt 

dl„,  _  4„  -  4 


5  COO 


1  (• 

J  * 

^  a=4^ 


*Equilibrium  phonon  intensity 
determined  by  Bose-Einstein  statistics 


dx 


X  A  0 


CO 


For  each  angle  (0),  solve  non-linear  equation  with  iterations  for 

—  Solving  for  Ioy 

-  Updating  Imo. 


n'COD 

0  / COS  OdcodCl 


Q=4  jt 


Internal  energy:  w(«  cT)  =  -^-^fi(of  D{(o)d(od€l=  j*  £ D —dcodQ 


Q=4  it 


Modeling  &  Solution  Procedure  using  COMSOL 


dl_ 

dt 


+  VCOS0 


dx 


-I 


T 


•  Use  a  built-in  feature  of  COMSOL,  “Coefficient  Form  PDEs”. 

•  The  spatial  domain  is  discretized  using  FE  mesh. 

•  The  angular  (momentum)  domain  is  discretized  using  Gaussian  quadrature 
points. 

•  For  each  angle  (0),  set  up  the  BTE  with  corresponding  coefficients  (ju^cosO) 
and  BCs  (Neumann  vs.  Dirichlet). 

•  Calculate  equilibrium  phonon  intensity  (/0)  by  numerical  integration  of  (/,)  using 
Gaussian  quadratures. 

•  Solve. 

—  Direct  solver  (UMFPACK). 

—  Max.  BDF  order  =  1 . 

•  Postprocess  and  visualize  the  results. 
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Details  of  Solution  Procedure 


Original  1-D  BTE:  (2L 

dt 


M  W  , 

VjLl - =  — -  ,  (jH-  COSO) 

dx  T„ 


t  x  yV 

•  Nondimensionalize  with  t  =  — ,  77  =  —  ,  Kn  =  — 


L 


L 


.  dI  r  &  r  r 

—r  +  Knju  —  +  I  =  L 


dt 


drj 


Split  into  (+)  and  (-)  directions: 

Discretize  angular  space  at 
Gaussian  quadrature  points: 


5/,+  r  m;  t+ 
dt  1  dr/ 


I-  =  I o  ,  (Mi>  0) 


dt  di 7 


iM,  <  0) 


/; 

_  &  rpA 

•> 

/: 

_<J  jiA 

O 

II 

0 

II 

1 

77=1  7T 

77=1 


After  FE  run,  postprocess: 


4(^)  = 


l 


q(t,rj)  =  2  n 


nSP /2  ”ff/2 

Z  w,7,+  +  Z  wf7r 

/=i  7=1 


”g p/2  /2 

7=1  7=1 


9 


Coefficient  Form  for  BTE 

(60  Finite  elements,  16  Gaussian  Points) 


File  Edit  Options  Draw  Physics  Mesh  Solve  Postprocessing  Multiphysics  Help 


D  &  U  #  V  l 


fet  A  A  #4  =  9*  H  Kl  Ci  #  $9 


BGeoml 


PDE,  Coefficient 
PDE,  Coefficient 
PDE,  Coefficient 
PDE,  Coefficient 
PDE,  Coefficient 
PDE,  Coefficient 
PDE,  Coefficient 
PDE,  Coefficient 
PDE,  Coefficient 
PDE,  Coefficient 
PDE,  Coefficient 
PDE,  Coefficient 
PDE,  Coefficient 
PDE,  Coefficient 
PDE,  Coefficient 


Form  (c2) 
Form (c3) 
Form (c4) 
Form (c5) 
Form (c6) 
Form (c7) 
Form (c8) 
Form (c9) 
Form  (dO) 
Form  (dl) 
Form  (d2) 
Form  (d3) 
Form (cl 4) 
Form  (d5) 
Form  (c!6) 


PDE,  Coefficient  Form  (c) 

Dependent  variables:  ul 
Default  element  type:  Lagrange  -  Quai 
Wave  extension:  Off 
Weak  constraints:  Off 


< 

1(0.901) 


A 

A 

2 

A 

A 

A 

A 

A 

A 

▲ 


t3 


Equation 

V'(-cVul  -  aul  4-  y)  4-  aul  4-  p'Vul  =  f 


Subdomains  Groups 


Subdomain  selection 


Group: 

[— ]  Select  by  group 
@  Active  in  this  domain 


Coefficients  init  Element  Weak 


PDE  coefficients 


Coefficient  Value/Expression  Description 

0  Diffusion  coefficient 


Io 


I 


|Kn*mui(l)~ 


Absorption  coefficient 
Source  term 
Mass  coefficient 
Damping/Mass  coefficient 
Conservative  flux  convection  coeff. 
Convection  coefficient 
Conservative  flux  source  term 


OK  Cancel  Apply  Help 


0  0.1  0.2  0.3  0.4  0.S  0.6  0.7  0.8  0.9 


Mesh  consists  of  15  elements. 
Mesh  consists  of  30  elements. 
Mesh  consists  of  60  elements. 


t~\ — i  r~~r 


Normal 


'Memory:  (116  / 122) 
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Nondim.  emissive  power,  eb 


Steady-State  Problem:  Analytic  vs.  Numerical  Solutions 


Emissive  power  -  Temperature 


Gradient 


J+ 

J  q\  J  ql 


Ae0=e0(T1  =  0)-e0(T]  =  l) 
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Nondimensional  heat  flux 


Steady-State  Problem:  Analytic  vs.  Numerical  Solutions 


Heat  flux 


Thermal  conductivity 


q 


q 


j+  -J 

J q\  J  q2 
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Nondimensional  temperature, 0 


Transient  Problem:  Effect  of  Spatial  Refinement 

(More  Finite  Elements) 


•  Refine  spatial  (x-)  direction  with  n  finite  elements  (n=15,  60,  120). 

•  Divide  angular  direction  with  16  Gaussian  points  (ngp=16). 


•  However,  it  does  not  solve  ray  effect. 


Nondimensional  temperature, 0 


Transient  Problem:  Effect  of  Angular  Refinement 

(More  Gaussian  Points) 

•  Refine  spatial  (x-)  direction  with  240  FE  elements. 

•  Divide  angular  direction  with  ngp  Gaussian  points  (ngp=4,8,16,32,64,128). 


•  Highly  refined  spatial  mesh  with  coarse 
angular  mesh  alleviates  solution. 
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Nondimensional  heat  flux,  q*  Nondimensional  temperature^ 


Results  of  BTE 

(Temperature  &  Heat  Flux  Distributions  with  Time  Increase) 


Nondimensional  coordinate,  r| 


Nondimensional  coordinate,  r) 


Nondimensional  coordinate,  r| 


Nondimensional  coordinate,  rj 


Nondimensional  coordinate,  rj 
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Nondimensional  coordinate,  rj 
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Nondimensional  heat  flux,  q*  Nondimensional  temperature^ 


Comparisons  of  FE,  HHTC,  BDE  vs.  BTE 

(Temperature  &  Heat  Flux  Distributions  for  3  Kn) 


Nondimensional  coordinate,  rj 


Nondimensional  coordinate,  rj 


Nondimensional  coordinate,  rj 


Nondimensional  coordinate,  rj 


0  0.2  0.4  0.6  0.8  1 

Nondimensional  coordinate,  rj  1  g 


Summary  and  Conclusions 


•  Nanoscale  simulation  is  conducted  for  phonon  heat  transfer  using 
Boltzmann  transport  equation. 

—  Phonons  for  dielectric,  thermoelectric,  semiconductor  materials. 

—  Electrons  for  metals. 

—  Gas  molecules  for  rarefied  gas  states. 


•  Numerical  solution  of  BTE  has  been  obtained  for  1-D  problem  (both 
steady-state  and  transient  problems). 


•  Temperature  and  heat  flux  distributions  from  the  nanoscale  simulation 
yield  completely  different  results  from  the  solutions  from  Fourier  and 
Cattaneo  equations. 

—  Thermal  conductivity  will  be  different,  too. 
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